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Abstract 

We present an exact mathematical transformation which converts a wide class of advection- 
diffusion equations into a form allowing simple and direct spatial discretization in all dimensions, 
and thus the construction of accurate and more efficient numerical algorithms. These discretized 
forms can also be viewed as master equations which provides an alternative mesoscopic interpre- 
tation of advection-diffusion processes in terms of diffusion with spatially varying hopping rates. 

PACS numbers: 02.70.Bf, 47.27.-i 
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I. INTRODUCTION 



Advection-diffusion equations (ADEs) describe a broad class of processes in the natural 
sciences. As their name implies, they provide a continuum (macroscopic) representation of 
systems whose underlying dynamics combines Brownian motion (diffusion) with some form 
of deterministic drift (advection). In this paper we shall consider ADEs of the general form 

dtp = V ■ DVp - V ■ pv . (1) 

The field p typically describes the number density of "particles" which, depending on the ap- 
plication, can range from electrons in a plasma, to chemical molecules advected in solution, 
to colloidal particles, to biological cells moving along chemical gradients. In principle the 
diffusion coefficient D(x, t) and the velocity field v(x, t) can depend on the density field p. 
An idea of the ubiquity of ADEs can be gauged from their diverse applications to traditional 
physics, soft matter systems, and biology. A small subset of examples are magnetic fusion 
plasmas Q], cosmic ray streaming 0, |J, electrons in weakly ionized gases ^, microemul- 
sions under shear fiow P|, chemical kinetics in driven systems 0, Q|, hydrodynamics and 
chemotaxis of bacterial colonies 0,0], phase field dynamics in directional solidification jlol |. 
and a wide array of tracer diffusion problems (for example 

It is generally not possible to analytically solve ADEs, especially since they often appear 
within sets of non-linear coupled equations. For this reason, great emphasis has been placed 
on numerical integration methods, typically based on finite differences. It has been found 
;hat the advection term, despite its apparent simplicity, is extremely troublesome to handle 



12| . There are two major challenges: stability, which can be improved using a range of 
implicit methods, and accuracy, which is a delicate issue, requiring the "best possible" form 
of spatial discretization. Regarding the issue of stability, many schemes are in use, such 
as the Crank-Nicholson/ ADI and fractional step methods Q], and the Lax-Wendroff 
method The issue of accuracy has received somewhat less attention with two spatial 



discretization schemes (and their immediate variants) commonly in use: these are the simple 
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171 ]. One of the main results of this 



Taylor expansion and the "upwind" scheme 

paper is the derivation of a new discretization scheme which is physically appealing, simple 
to apply in all dimensions, and more accurate than those currently in use. 
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II. A SIMPLE EXAMPLE 



To begin, let us present the key idea in the context of a simple ADE, namely, a one- 
dimensional system with a velocity function proportional to the spatial derivative of a scalar 
potential (j){x,t). Thus, we consider 



dtp = Dodlp - ad^{pd^(f)) , 



(2) 



where Dq and a are constants. 

Most numerical algorithms designed to integrate an equation such as (j2I) treat the diffusion 



and advection terms separately 
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17| . The difficulties arise in finding a discretization 



for the latter term. In doing so, two fundamental properties of the equation must be exactly 
maintained. These are the non-negativity of p and its spatial conservation: 



J dx p{x, t) = const . 



(3) 



As an illustration, let us write down a common spatial discretization using simple Taylor 
expansion, which is used both for ex plic it Euler schemes and as the basis for more 



advanced implicit algorithms 
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isl |l6[. We replace the continuous function p{x,t) by a 



set of functions {piit)} defined on a regular grid with lattice spacing h. The equation of 
motion for pi is written using centered spatial derivatives: 
dpi Do , 



-(Pi+i +P: 



dt ^^^^^ ' "^^^ 4/^2 

It is noteworthy that this simple scheme requires knowledge of the scalar field at next- 
nearest neighbor grid points rather than neighboring grid points. For future reference we 
rewrite this discrete equation in the following manner: 



2pi 



a 



(pi) - Pi-l(0i - 0i_2)] 



(4) 



dpi 
dt 



Pi+i 



a 

I' 



+ Pi~l 



a 



H-2 



(5) 



which we shall hereafter refer to as the "linear centered discretization" (LCD) (and which 



resembles the backward Euler scheme used for simple advection problems 

We now turn to a new discretization scheme which emerges from a simple mathematical 
transformation of the ADE Q. Defining 7 = ajlD^ it can be verified by direct differentia- 
tion that Eq.Q may be written as 



dtp = Do [e""^ dlipe-""^) - e-""^p dl{e''*)] . 



(6) 



A similar transformation involving exponential functions is known for Fokker-Planck equa- 
tions The simple ADE given in Q can indeed be formally interpreted as such an 
equation, although the physical origin is quite different. We will shortly be considering 
more general ADEs in which the diffusion coefficient and velocity function can be functions 
of the density p. Clearly then the simple correspondence with Fokker-Planck equations 
breaks down, although we are still able to achieve a transformation of the kind given above. 
The crucial feature of Eq.® is that spatial derivatives only enter in the form of a second 
derivative dl which is straightforward to discretize. Using the simplest such discretization 
we immediately have 



There are a number of points to make concerning this equation. First, in contrast to 
the LCD (^, the scalar field appears in a non-linear fashion, and is sampled at nearest- 
neighbor positions. Second, the new equation is of the same form as a master equation 



|20|. Within this analogy one can think of pi as the probability that a fictitious 
particle is located at grid position i. The transition rate for the particle to hop from grid 
point i to a neighboring point j is of the Arrhenius form 

Wi^, = {Do/h') exp[-7(0, - 0,)] . (8) 

Given this formal analogy with a master equation for a probability function, one immedi- 
ately sees that Eq.Q exactly maintains conservation of the function p (normalization of 
probability) and its non-negativity. Due to this analogy we hereafter refer to Eq.((7I) as the 
"master equation discretization" (MED). 

Our numerical work (see section V) shows that the MED is more accurate than the LCD 
and other popular discretizations. To appreciate the underlying reason for this, it is helpful 
to consider the case of 7(50 -C 1 in which case we can expand the exponential functions in 
Eq.(I7j) to first order. One then finds 



dpi _ 1 
dt ~ 



Pi+i 



a 



+ Pi-i 



a 



2Do- 2(20.-0m-0.-i) \ ■ (9) 



Comparison of this form with Eq.(jS)) gives useful insight into the potential weakness of the 
LCD. Namely, it neglects an important curvature term in the scalar field. In fact, this 



omission is directly related to artificial (or "numerical" ) diffusion, which is a common failing 



of other discretization schemes, most notably, the "upwind" scheme |12l. Il6l Il7|. The linear 
scheme given above in Eq.Q can of course be regarded as one of many possible linear 
discretizations, but without the derivation given here one would have no a priori reason to 
prefer it over forms such as the LCD, since they both have non-vanishing second-order errors 
in space. Continuing the expansion of the exponential terms in powers of a yields crucial 
non-linear corrections to Eq. Q which have no analogy within linear discretization schemes. 
As shall be seen below, the MED is easily formulated for the d- dimensional extension of 
Eq. (j2I) as well as for a range of more general ADEs. 



III. THE GENERAL CASE 



Consider the general ADE in (i-dimensions given in Eq.(0). We shall now proceed to 
transform this equation into a form amenable to the MED. In one dimension we shall find 
that this is possible for general functions D and v. In higher dimensions the vectorial nature 
of the velocity field will place a constraint on the transformation. 

Let us introduce two scalar functions /(x, t) and g{'x,t) defined via the relations 

D = fg, (10) 
V = ^7V/ - fVg . (11) 

Then the ADE has the explicit form 

= V • [fgVp] - V ■ [pigVf - fVg)] . (12) 

By direct differentiation one can show that this equation may be rewritten as 

dtp = fV'igp)-gpV'f. (13) 

Once again, we see that the spatial derivatives appear only as Laplacians, which allows us 
to immediately write down a simple discrete form. Let us define the discrete Laplacian via 

where the sum is over nearest neighbors j of the grid point i, which corresponds to the 
continuum position x. Then the MED corresponding to Eq. (jl3p is 

dtpi = Wj^i Pi - ^i^j pi] ' (15) 
j 
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where the transition rate for "hopping" from site i to site j is 

Wi^, = f, g,/h^ . (16) 

Having formulated the MED in this general manner, let us examine some particular cases. 
We stress that once the functions / and g are determined the discrete algorithm is completely 
defined via the transition rate given above. 

First, we consider one dimension. In this case it is possible to integrate Eqs. and 
(11111 exactly to find the necessary auxiliary functions / and g in terms of the physically 
relevant diffusion coefficient and velocity. One finds 

X 

f{x,t) = Cv/^(M)exp(5), S{x,t) = \j dx' , (17) 

with g then given trivially from pOj) . The transition rate is easily evaluated from ()16|1 to 
give 

W^.-., = ^^^ exp[-(5,-5,)] . (18) 

A non-trivial application of this general solution would be advection-diffusion in the kinetic 
theory of gases where the diffusion coefficient is non-constant, and actually depends on 
the density as -D oc 1/p j2l|. In higher dimensions a general solution for / and g is not 
possible. Solvable cases will rely on special conditions for D and v reminiscent of the 
potential conditions for the existence of steady-state solutions to the multi-variate Fokker- 
Planck equation 2ol |. 

For many problems the diffusion coefficient is constant (-Do) and the velocity function 
is associated with a scalar potential via v = q;V0. In these cases, the analysis leading to 
Eq.(jHI) is easily generalized to d dimensions and one finds the discrete equation (|T5|l with 



Wi^, = (Do/h^) exp[-7(0, - 0,)] , (19) 

where we remind the reader that 7 = a/2DQ. As found in one dimension, this scheme 
includes important curvature terms, even within a linear approximation, which are absent 
in conventional LCD algorithms. Numerical analysis shows such terms to be essential in 
regions where has maxima or minima. 

The MED scheme encapsulated in Eqs. ()15|) and p9p can be used to model more com- 
plicated ADEs in which there is non-linear feedback. An interesting example of this is the 
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continuum theory of group dynamics, in w 
nism is imposed via the velocity potential 



lich a non-linear and non-local feedback mecha- 



23|. In particular one has 



(20) 



where V is analogous to a potential, and is responsible for long-range attraction and short- 
range repulsion of individuals. If is a Dirac (5-function then oc p. Such models are used 
to describe density-dependent dispersal in population dynamics j22] and have recently been 
shown to arise from excluded volume effects in models of interacting cellular systems 241. 



A second well-known example is the Keller-Segel model for chemotactic motion j2|. Here, 
the potential represents the chemoattractant concentration field and is coupled to the cell 
density field p via 

dt4> = z/VV -X^ + pp , (21) 

where z/, A and f3 are the diffusion constant for the chemical field and its rate of degra- 
dation and production respectively. This equation is easily discretized and the resulting 
discrete chemical concentration field may be inserted into the transition rate (fTI?|) allowing 
a straightforward scheme for integration of the cell density. 



IV. FINE-TUNING THE MED ALGORITHM 



From numerical investigations (see next section) we have found that the MED is generally 
far more accurate than both the LCD and upwind schemes. In regions where the velocity 
function has strong spatial variation, the MED does an excellent job in predicting the correct 
density even for grid scales approaching the scale of variation of the velocity. However, in 
the "simpler case" when dynamics are dominated by advection in a region of quasi- const ant 
velocity, the MED fares less well. This problem can be traced back to the exponential 
weights yielding, in regions of constant velocity, an over-estimated drift velocity. In terms 
of a hopping process, the bias in hopping rates between neighboring sites is proportional to 
sinh(750), whereas the correct drift velocity is simply proportional to 7(50. 

We discuss here two straightforward extensions to MED which alleviate this problem, but 
also lead to slightly less accurate algorithms in the "non-trivial" regions where the velocity 
is strongly varying. Both extensions amount to a renormalization of the hopping rates. 
An ideal algorithm would be a hybrid, using the original MED and either of the following 
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extensions in appropriate regions. We will not discuss such hybrid schemes here since their 
form will be strongly dependent on actual applications. 

For simplicity let us consider again the one- dimensional ADE given in Eq.Q. The MED 
scheme for this case in given in Eq.((7j), where the transition rate from site i to neighboring 
site j has the explicit form 

W^^j = (Do/h^) exp[-7(0, - 0,)] . (22) 

It is clear from ()22|) that the effective drift velocity arising from the bias in hopping rates 
between i and j is 

Ves = h{Wi^j - Wj^.i) = (2Do//i)sinh[a((/), - 0.)/2Do] , (23) 

where we have reinstated 7 = a/2DQ for clarity. The correct drift velocity between these 
two points is simply a{(j)j — <f)i)/h which is recovered if the grid scale is small (or else the 
velocity potential is slowly varying). 

In order to correct the MED algorithm one may either renormalize the effective diffusion 
coefficient (which is the pre-factor of the exponential weight) or else renormalize the param- 
eter 7 which appears in the argument of the exponential. In the former case one has, on 
fitting the drift velocity to its correct value, the effective diffusion coefficient 

- ^°sinh(a50/2Do) ^ ^ 

which leads to the MED transition weight taking the "Fermi-Dirac" (FD) form 

W- = «(0. - 0,)/Do 

UV exp[a(0i - 0,)/Z^o] - 1 ■ ^ ^ 

The alternative is to correct the drift velocity by adjusting 7, which leads to 

7., = A-nh- ( . (26) 



6^ \2Do/ 

Writing the inverse hyperbolic function in terms of a logarithm leads to the MED transition 
rate taking the "square root" (SR) form 



2Dn 



1/2 



ft(0i - 0, 
2Dn 



(27) 



Numerically one finds that the FD form ()25|) is generally more accurate than the SR form 
and that both are superior to the LCD and upwind schemes. As already mentioned, 
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the original MED scheme defined by Eq.()22|) is the best of all the schemes described when 
the velocity field is strongly varying, and/or during asymptotic relaxation of the density 
field to its steady-state. 



V. NUMERICAL WORK 

We have made a careful numerical analysis of the simple one-dimensional ADE given in 
Eq.Q, along with its two-dimensional extension. Since we wish to gauge the accuracy of 
our new scheme, we have compared the MED scheme ([7j), and its variants (the MED(FD) 
given in the MED(SR) given in (gZj), and the linearized MED, denoted by MED(LIN), 
given in (jH))), with both the LCD and upwind schemes In one dimension we use a static 
velocity potential given by <f){x) = [1 + cos{2Tmx/L)]/2 with n = 16 and L = 12.8. The 
initial density function is taken to be uniform in the region x G (—3, 3) and zero otherwise. 
The density is normalized to unity and periodic boundary conditions are enforced. This set- 
up provides a challenging test of all the schemes since the velocity field is a strongly varying 
function of position. Furthermore, we challenge the methods by using the parameter values 
-Do = 1-0 and a = 5.0 (Figiire 1) and a = 20.0 (Figure 2), which correspond to moderate to 
high grid Peclet numbers p| at the grid scales of interest. Here, the largest Peclet numbers 
are approximately given by 2ah and so vary between 0.25 and 8 for the data shown in 
Figures 1 and 2. The dynamics consists of a rapid transient phase where the density field 
adapts to the periodic structure of the velocity field, followed by a slower relaxation towards 
the steady-state. Thus, the numerical analysis probes each scheme's ability to track rapid 
advective motion and diffusive relaxation around maxima and minima of the velocity field. 

In order to assess the accuracy of the methods we first run all schemes at a very small grid 
size oih = 0.00625, using an explicit temporal scheme with 5t = 10~^. Very good agreement 
is found among all the schemes and the solution is denoted "exact." We then run all the 
schemes at larger grid scales using 6t = 10~'^, and dynamically compare the approximate 
solutions with the exact one. This is gauged using the relative error, which is defined via 

Y^[Pi(t) - Pi,cxact(t)]^ 

^^^^ = (tV- ^^^^ 

/ J exact V''/ 
i 

Note, that St is chosen small enough such that any differences between our first-order tem- 
poral discretization for LCD and second-order schemes (in the temporal dimension) such as 



Crank-Nicholson or Lax-Wendroff are negligible. Figures l(a)-(e) show E{t) for grid scales 
h = 0.025, 0.05, 0.1, 0.2, and 0.4 respectively, for a = 5.0. The entire dynamical evolution 
up to the steady state is shown. In the first four panels we clearly see that the MED and 
its (nonlinear) variants give a relative error approximately 10 times less than the LCD and 
UW schemes. (Note UW does not appear in 1(a) since its error is too large to be usefully 
included in the figure.) The relative errors of all the schemes increases roughly by a factor 
of 10 as the grid scale is doubled. Panel 1(e) shows the breakdown of all the schemes at the 
scale h = 0.4 which is comparable to the period of the velocity field. By "breakdown" we 
mean a relative error of 10% or more. To give an idea of the spatial form of the density field 
near the steady-state we show in Figure 1(f) the exact density profile in a peripheral region, 
along with the LCD and MED (FD) at a grid scale of h = 0.2 for comparison. Note the 
LCD fails to capture the magnitude of the maximum density, and also becomes negative at 
some grid points. 

In a similar fashion, figures 2(a)-(d) show E{t) for h = 0.025, 0.05, 0.1, and 0.2 respec- 
tively, for a = 20.0. As before the non-linear MED schemes perform far better than the 
LCD and UW, meaning the relative error is roughly 10 times smaller for a given grid scale. 
Note also that the MED(FD) and MED(SR) algorithms perform better than MED during 
the transient period, as expected. All schemes break down for h = 0.2. In Figure 2(e) we 
show the exact density profile close to the steady-state, compared with the MED and LCD 
schemes for h = 0.1. Again, the LCD shows negative values and fails in the vicinity of the 
density peaks. Figure 2(f) is the same except the UW scheme is compared to the MED. The 
UW scheme is designed to give non-negative densities, but has high (artificial) "numerical 
diffusion" which inflate the width of the density peaks. 

We have performed an exactly analogous numerical examination in two dimensions. We 
integrated the 2D generalization of Eq. using the potential 0(x, ?/) = [l+cos(27rnx/L)][H- 
cos(27rny/L)]/4 with n = 16 and L = 12.8. We take Dq = 1.0 and a = 10.0. The initial 
density function is uniform in a disk of radius 3.0 and zero otherwise, and again normalized 
to unity. The "exact" density profile is evaluated using h = 0.0125 and 6t = 0.25 x 10~^. 
The two-dimensional extensions of all six schemes are integrated for grid scales of h = 
0.025, 0.05, 0.1, and 0.2 using 6t = 10~'^. The relative error E{t) for these cases is shown in 
Figure 3(a)-(d), for a time period encompassing the initial rapid adaptation to the potential 
followed by the early stages of relaxation to the steady-state. As with one dimension, the 
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MED and its (non-linear) variants perform far better than the LCD and UW, with the pure 
MED scheme performing best at later times. All schemes break down for h = 0.2. Direct 
comparison of the exact density profile, MED, and LCD is given in Figures 3(e) and (f), 
for h = 0.05 and h = 0.1 respectively, along a one-dimensional cut {y = 0) in a peripheral 
region of the density. The MED shows excellent agreement, especially in the vicinity of the 
density peaks. The LCD fails in the vicinity of the density peaks as expected. 

From this and similar numerical work we have concluded that the MED and its (non- 
linear) extensions are superior spatial discretization schemes compared to the LCD and 
upwind schemes. The MED works especially well in regions of large variation in the velocity 
potential. Generally speaking, for a given error tolerance, the MED and variants allow one 
to use grid scales at least two times larger than traditional schemes, which translates into 
a saving of at least a factor of 4 and 8 in computational cost for two and three dimensional 
numerical analyses. 

VI. DISCUSSION AND CONCLUSIONS 

We end with some remarks on the non-linear transition rates of the MED. In most appli- 
cations the ADEs represent processes for which there is no underlying lattice (e.g. cosmic 
ray diffusion 3] or chemotactically moving cells j^). When one discretizes the continuum 
ADE one must therefore not regard the lattice version as "more fundamental" or "more mi- 
croscopic." It is simply a mathematical analog of the original equation and identical in the 
limit of the lattice spacing being taken to zero. This is a different situation to that found for 
many models arising from solid state physics in which there is an underlying crystal lattice, 
and for which the discrete equation can often be regarded as more fundamental (or, at least, 
more microscopic) than continuum models. Although the hopping process encapsulated by 
the MED cannot be viewed as the underlying microscopic dynamics, it is interesting that 
ADEs can be accurately modeled by a process in which diffusion and advection are non- 
linearly combined in Arrhenius transition rates. Figure 4 summarizes our understanding of 
the algorithmic connections between ADE and the MED discretization, in which a given 
ADE typically arises from a mean-field approximation of a microscopic stochastic process 
which is not constrained by a lattice. 

Pragmatically one wishes to impose a "large" lattice scale for numerical efficiency, while 
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avoiding the loss of accuracy. Algorithms which remain accurate for larger lattice scales yield 
great computational speed-up in higher dimensions, since the number of required grid points 
(and hence computer operations) scales as h~'^. We find that our new scheme typically allows 
grid scales between 2-4 times larger than traditional schemes, which in three dimensions 
allows a potential speed-up in computation of one or two orders of magnitude. Naturally, 
our improved spatial discretizations can be used in more advanced algorithms which use 
implicit temporal methods and/or adaptive spatial grids. 

In conclusion we have shown that a wide class of advection-diffusion equations can be ex- 
actly rewritten in a form which immediately allows a direct and simple spatial discretization 
in all dimensions. Our new discrete forms contain important non-linear terms, which when 
linearized are seen to be related to the curvature of the velocity potential, such terms being 
absent in commonly used discretization schemes. We have shown explicitly that these curva- 
ture effects are essential for accurate integration of ADEs, both in one and two dimensions, 
and allow simple algorithms to be constructed which are accurate for grid scales up to the 
size of spatial variation in the velocity field. We estimate that our new algorithm may allow 
a speed-up of ADE computation by factors of 10 or more in three dimensions due to the 
increased grid scale one can impose. The fact that ADE can be recast as master equations 
also yields interesting physical insight into their dynamics - namely that at mesoscopic scales 
the processes of diffusion and advection may be modeled as a non-linear combination within 
Arrhenius-like transition rates. 

The authors gratefully acknowledge partial support from NSF award DEB-0328267. 
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Figure Captions 

Figure 1: Data from numerical integration of Eq.Q using various schemes in one dimension, 
with Dq = 1.0 and a = 5.0. The particular form of the velocity potential and the initial 
density profile are described in section V. The time step is 6t = 10~^. Figures 1(a), (b), 
(c), (d), and (e) show the relative error as a function of time for grid scales oi h = 
0.025, 0.05, 0.1, 0.2, and 0.4 respectively. The methods used are upwind (UW), LCD 
dH), linearized MED ©, MED (0), 'Fermi-Dirac" version of MED and "square-root" 
version of MED ()27|). Figure 1(f) compares the exact density profile in the peripheral region 
X G (2,3.6) with both the MED(FD) scheme and the LCD scheme at time t = 0.1 using 
h = 0.2. In Figures 1-3, time is measured in units of 6t, space in units of h, and the density 
in dimensionless units. 

Figure 2: Same as Figure 1, but with a = 20.0. Figures 2(a), (b), (c), and (d) show 
the relative error ()28|) as a function of time for grid scales of h = 0.025, 0.05, 0.1, and 
0.2 respectively. Figure 2(e) compares the exact density profile in the peripheral region 
X G (2, 3.6) with both the MED scheme and the LCD scheme at time t = 0.02 using h = 0.1. 
Figure 2(f) is the same as 2(e) but compares the exact profile with both MED and UW. 

Figure 3: Data from numerical integration of the two-dimensional generalization of Eq.(j21) 
using various schemes, with Dq = 1.0 and a = 10.0. The particular form of the velocity 
potential and the initial density profile are described in section V. The time step is 6t = 10~^. 
Figures 1(a), (b), (c), and (d) show the relative error (j^Hj) as a function of time for grid 
scales oi h = 0.025, 0.05, 0.1, and 0.2 respectively. The methods used are two dimensional 
generalizations of upwind (UW), LCD (jH), linearized MED ©, MED ©, "Fermi-Dirac" 
version of MED ()25|). and "square-root" version of MED ()27p . Figure 3(e) compares the 
exact density profile along a cut {y = 0) in the peripheral region x G (2, 3.6) with both the 
MED scheme and the LCD scheme at time t = 0.01 using h = 0.05. Figure 3(f) is the same 
as 3(e) except that a larger grid scale of /i = 0.1 is used. 

Figure 4: A schematic diagram summarizing the relationships between various descrip- 
tions of advection-diffusion processes. The MED is a useful mesoscopic description in terms 
of Arrhenius hopping rates, rather than a reflection of the underlying dynamics. 
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